rm(list = ls(all.names = TRUE))
gc()
set.seed(211917)
options(scipen=999)

packages <-c("tidyverse","foreign","stargazer","naniar","ri2")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR WORKING DIRECTORY HERE")

# load data
dat <- read.dta("./data.dta")

#---------------#
# Data cleaning #
#---------------#

dat$q6_years_married[dat$q6_years_married == -9] <- 0

## select all possible covariates and outcome at first
dat <- dat %>% 
  select(q30_combined, treatment, holy_endorser, secular_endorser, respondent_muslim, respondent_male, enumerator_muslim, enumerator_male, q1_age, q2_male, q3_kamba, q3_kikuyu, q3_luo, q3_somali, q3_other, q4_religion_christian, q4_religion_muslim, q5_convert, q6_years_married, q6_married, q7_has_job, q8_income_last_month, q8_has_income, q9_chance_become_rich, q10_plans_to_vote, q11_govt_represents_interest, q12_feels_on_loosing_side, q13_a_knows_emigrant, q13_b_knows_emigrant_som_eri, q15_witnesses_interrel_violence, q16_witnesses_viol_musl_govt, q17_friend_family_died, q18_lost_job, q18_arrested, q18_house_of_worship_raided, q18_stopped_speaking_parents, q18_love_problems, q18_friend_left_country, q18_sum, q19_rlts_mother, q19_not_raised_by_mom, q19_mother_dead, q20_rlts_father, q20_not_raised_by_father, q20_father_dead, q21_respect_friends_family, q22_gender, q22_religion, q22_tribal, q22_national, q22_youth, q23_religious_attendance, q23_frequents_house_worship, q24_deontology, q25_payment_acceptable, q26_payment_ethical, q27_radicalization, q33_distance)

## check for missings
vis_miss(dat)

## following variables seems problematic, i.e., shouldn't impute at this level
table(is.na(dat$q13_a_knows_emigrant))
table(is.na(dat$q19_rlts_mother))
table(is.na(dat$q20_rlts_father))
table(is.na(dat$q23_religious_attendance))

## remove problematic vars
dat <- dat %>% select(-q13_a_knows_emigrant, -q19_rlts_mother, -q20_rlts_father, -q23_religious_attendance)

## check missing again:
vis_miss(dat) 

## remove one missing treatment indicator (can't be imputed)
dat <- dat[which(dat$holy_endorser>-1),]

## impute rest of missings with mean 
for(i in 1:ncol(dat)){
  dat[is.na(dat[,i]), i] <- mean(dat[,i], na.rm = TRUE)
}
rm(i)

## 
dat <- dat[complete.cases(dat), ]
dat$q30_combined_s <- scale(dat$q30_combined)

christian_d <- dat[which(dat$respondent_muslim==0),]
muslim_d <- dat[which(dat$respondent_muslim==1), ]
killing_d <- dat[which(dat$treatment==1),]
suicide_d <- dat[which(dat$treatment==2),]

model1 <- lm(q30_combined ~ holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=dat)
model2 <- lm(q30_combined ~ holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=killing_d)
model3 <- lm(q30_combined ~ holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=suicide_d)
model4 <- lm(q30_combined ~ holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=christian_d)
model5 <- lm(q30_combined ~ holy_endorser + secular_endorser + respondent_male + q1_age + respondent_muslim + q5_convert + q6_married + q6_years_married + q7_has_job + q8_has_income + q8_income_last_month + q3_kamba + q3_kikuyu + q3_luo + q3_somali + enumerator_muslim + enumerator_male, data=muslim_d)

vars <- c("Religious Endorser", "Secular Endorser",	"Male",	"Age","Muslim",	"Convert",	
          "Married", "Years Married", "Employed","Any Income", "Income (KES)",	
          "Tribe Kamba", "Tribe Kikuyu", "Tribe Luo", "Tribe Somali",	
          "Enumerator Muslim", "Enumerator Male")

stargazer(model1, model2, model3, model4, model5,
          covariate.labels = vars,
          keep.stat = c("n"),
          digits = 3,
          dep.var.caption = "Support for Violence",
          dep.var.labels = c("Full","No Killing","No Suicide","Christian","Muslim"),
          column.labels = c("Full sample","No Killing","No Suicide"),
          label = "tab:main",
          title = "Treatment effects on violent attitudes",
          out = "PUT YOUR FILE PATH HERE")





